"""
GMSH操作模块
包含所有与GMSH相关的网格处理函数
"""

import os
import json
import gmsh
import math
import hashlib


def create_feature_id_from_gmsh(center, area):
    """
    从Gmsh面数据创建特征标识符（与OCC保持一致）
    """
    # 创建特征字符串（与OCC格式保持一致）
    feature_str = f"{area:.6f}_{center[0]:.6f}_{center[1]:.6f}_{center[2]:.6f}"
    return hashlib.md5(feature_str.encode()).hexdigest()[:16]


def find_matching_face_with_features(gmsh_center, gmsh_area, feature_mapping, center_tol=1e-3, area_tol=0.1):
    """
    在特征映射中查找匹配的面（考虑容差）
    """
    target_feature_id = create_feature_id_from_gmsh(gmsh_center, gmsh_area)
    
    # 精确匹配
    if target_feature_id in feature_mapping:
        return feature_mapping[target_feature_id]
    
    # 容差匹配（如果精确匹配失败）
    for feature_id, data in feature_mapping.items():
        stored_center = data['center']
        stored_area = data['area']
        
        # 计算中心距离
        center_distance = math.sqrt(
            (gmsh_center[0] - stored_center[0])**2 +
            (gmsh_center[1] - stored_center[1])**2 + 
            (gmsh_center[2] - stored_center[2])**2
        )
        
        # 计算面积相对误差
        area_error = abs(gmsh_area - stored_area) / max(gmsh_area, stored_area)
        
        if center_distance < center_tol and area_error < area_tol:
            print(f"容差匹配: 距离={center_distance:.6f}, 面积误差={area_error:.6f}")
            return data
    
    return None


def load_feature_and_meshsize_data(faceid_to_features_file, faceid_to_meshsize_file):
    """
    加载faceid到features和meshsize的映射数据
    """
    import os
    
    # 加载faceid到features的映射
    with open(faceid_to_features_file, 'r', encoding='utf-8') as f:
        faceid_to_features = json.load(f)
    
    # 加载faceid到meshsize的映射
    with open(faceid_to_meshsize_file, 'r', encoding='utf-8') as f:
        faceid_to_meshsize = json.load(f)
    
    print(f"加载了 {len(faceid_to_features)} 个faceid的特征数据")
    print(f"加载了 {len(faceid_to_meshsize)} 个faceid的分网配置")
    
    return faceid_to_features, faceid_to_meshsize


def setup_gmsh_environment():
    """
    初始化GMSH环境
    """
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)
    

def configure_mesh_parameters():
    """
    配置网格生成参数
    """
    # 设置网格生成参数，确保几何不变形
    gmsh.option.setNumber("Mesh.Algorithm", 6)  # 使用Frontal-Delaunay算法
    gmsh.option.setNumber("Mesh.Algorithm3D", 1)  # 使用Delaunay算法
    gmsh.option.setNumber("Mesh.RecombineAll", 1)  # 重新组合所有网格
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)  # 使用复合算法


def match_face_with_features(face_entity, faceid_to_features, faceid_to_meshsize, face_index):
    """
    为单个面进行几何特征匹配
    """
    dim, tag = face_entity
    
    try:
        # 获取面的几何特征
        center = gmsh.model.occ.getCenterOfMass(dim, tag)
        area = gmsh.model.occ.getMass(dim, tag)
        
        # 在faceid-to-features数据中查找匹配的面
        matched_data = find_matching_face_with_features(center, area, faceid_to_features)
        
        if matched_data:
            # 找到匹配，使用对应的faceid
            matched_faceid = None
            for faceid, features in faceid_to_features.items():
                if features.get('center') == matched_data.get('center') and features.get('area') == matched_data.get('area'):
                    matched_faceid = faceid
                    break
            
            if matched_faceid:
                # 根据faceid获取分网尺寸
                mesh_size = faceid_to_meshsize.get(matched_faceid, 1.0)
                
                # 设置分网尺寸
                gmsh.model.mesh.setSize([face_entity], mesh_size)
                
                # 创建物理组（使用匹配的faceid作为标识）
                physical_group = gmsh.model.addPhysicalGroup(dim, [tag])
                gmsh.model.setPhysicalName(dim, physical_group, matched_faceid)
                
                print(f"面 {tag}: 几何特征匹配成功, faceid={matched_faceid}, 分网尺寸={mesh_size}")
                return True, matched_faceid
            else:
                # 未找到对应的faceid，使用默认配置
                default_size = 1.0
                gmsh.model.mesh.setSize([face_entity], default_size)
                physical_group = gmsh.model.addPhysicalGroup(dim, [tag])
                gmsh.model.setPhysicalName(dim, physical_group, f"unmatched_{face_index}")
                print(f"面 {tag}: 几何特征匹配但未找到faceid，使用默认尺寸 {default_size}")
                return False, f"unmatched_{face_index}"
        else:
            # 未找到匹配，使用默认配置
            default_size = 1.0
            gmsh.model.mesh.setSize([face_entity], default_size)
            physical_group = gmsh.model.addPhysicalGroup(dim, [tag])
            gmsh.model.setPhysicalName(dim, physical_group, f"unmatched_{face_index}")
            print(f"面 {tag}: 未找到几何特征匹配，使用默认尺寸 {default_size}")
            return False, f"unmatched_{face_index}"
            
    except Exception as e:
        print(f"处理面 {tag} 时出错: {e}")
        # 出错时使用默认配置
        default_size = 1.0
        gmsh.model.mesh.setSize([face_entity], default_size)
        physical_group = gmsh.model.addPhysicalGroup(dim, [tag])
        gmsh.model.setPhysicalName(dim, physical_group, f"error_{face_index}")
        return False, f"error_{face_index}"


def perform_feature_matching(iges_file_path, faceid_to_features, faceid_to_meshsize):
    """
    执行几何特征匹配流程
    """
    # 加载IGES文件
    gmsh.open(iges_file_path)
    
    # 获取所有面
    face_entities = gmsh.model.getEntities(2)
    
    print(f"步骤1: 为 {len(face_entities)} 个面进行几何特征匹配...")
    
    matched_count = 0
    unmatched_count = 0
    
    # 为每个面进行几何特征匹配
    for face_index, entity in enumerate(face_entities):
        is_matched, faceid = match_face_with_features(entity, faceid_to_features, faceid_to_meshsize, face_index)
        if is_matched:
            matched_count += 1
        else:
            unmatched_count += 1
    
    print(f"几何特征匹配统计: {matched_count} 个面匹配成功, {unmatched_count} 个面未匹配")
    
    return face_entities


def generate_and_save_mesh(output_dir):
    """
    生成网格并保存
    """
    import os
    
    # 生成网格（2D网格）
    gmsh.model.mesh.generate(2)
    
    # 保存网格文件
    output_mesh_file = os.path.join(output_dir, "gmsh_meshed_model.msh")
    gmsh.write(output_mesh_file)
    print(f"Gmsh网格已保存到: {output_mesh_file}")
    
    return output_mesh_file


def convert_mesh_to_vtu_with_meshio(mesh_file_path, output_dir):
    """
    使用meshio库将GMSH网格转换为VTU格式并添加FaceID
    """
    import os
    import meshio
    import numpy as np
    
    print("步骤2: 使用meshio库转换GMSH网格为VTU格式...")
    gmsh_vtu_file = os.path.join(output_dir, "gmsh_model_with_face_ids.vtu")
    
    # 使用meshio读取.msh文件
    mesh = meshio.read(mesh_file_path)
    print(f"已读取GMSH网格文件: {mesh_file_path}")
    
    # 获取所有物理组
    physical_groups = gmsh.model.getPhysicalGroups()
    
    # 创建物理组到FaceID的映射
    physical_group_to_faceid = {}
    for dim, tag in physical_groups:
        name = gmsh.model.getPhysicalName(dim, tag)
        physical_group_to_faceid[tag] = name
    
    # 为每个单元分配FaceID
    face_ids = []
    
    # 获取单元到物理组的映射
    element_physical_groups = {}
    for dim, tag in physical_groups:
        # 获取物理组包含的实体
        entities = gmsh.model.getEntitiesForPhysicalGroup(dim, tag)
        for entity in entities:
            # 获取实体包含的单元
            element_types, element_tags, _ = gmsh.model.mesh.getElements(dim, entity)
            for elem_type, elem_tags in zip(element_types, element_tags):
                for elem_tag in elem_tags:
                    element_physical_groups[elem_tag] = physical_group_to_faceid[tag]
    
    # 为每个单元分配FaceID
    for cell_block in mesh.cells:
        for idx in range(len(cell_block.data)):
            # 默认FaceID为0
            face_id = 0
            
            # 尝试从物理组映射中获取FaceID
            if idx < len(element_physical_groups):
                face_id = list(element_physical_groups.values())[idx % len(element_physical_groups)]
            
            face_ids.append(face_id)
    
    print(f"已为 {len(face_ids)} 个单元分配FaceID")
    
    # 添加FaceID作为单元数据
    face_ids_array = np.array(face_ids, dtype=np.int32)
    mesh.cell_data["FaceID"] = [face_ids_array]
    
    # 清除cell_sets以避免meshio转换错误
    if hasattr(mesh, 'cell_sets'):
        mesh.cell_sets = {}
    
    # 写入VTU格式
    meshio.write(gmsh_vtu_file, mesh)
    print(f"GMSH VTU模型已保存到: {gmsh_vtu_file}")
    
    return gmsh_vtu_file


def convert_mesh_to_vtk_fallback(mesh_file_path, output_dir):
    """
    备用方法：使用VTK将GMSH网格转换为VTK格式
    """
    import os
    import vtk
    import numpy as np
    
    print("警告: meshio库未安装，使用备用VTK转换方法")
    gmsh_vtk_file = os.path.join(output_dir, "gmsh_model_with_face_ids.vtk")
    
    # 获取网格数据
    node_tags, node_coords, _ = gmsh.model.mesh.getNodes()
    element_types, element_tags, element_node_tags = gmsh.model.mesh.getElements()
    
    # 创建VTK数据结构
    points = vtk.vtkPoints()
    polys = vtk.vtkCellArray()
    face_id_array = vtk.vtkIntArray()
    face_id_array.SetName("FaceID")
    
    # 添加节点
    coords = np.array(node_coords).reshape(-1, 3)
    for i in range(coords.shape[0]):
        points.InsertNextPoint(coords[i, 0], coords[i, 1], coords[i, 2])
    
    # 创建节点标签到索引的映射
    node_tag_to_index = {tag: idx for idx, tag in enumerate(node_tags)}
    
    # 处理所有单元类型
    cell_count = 0
    for i in range(len(element_types)):
        elem_type = element_types[i]
        elem_tags = element_tags[i]
        elem_node_tags = element_node_tags[i]
        
        # 检查数据格式
        if hasattr(elem_tags, '__iter__') and hasattr(elem_node_tags, '__iter__'):
            # 处理每个单元
            for j in range(len(elem_tags)):
                elem_tag = elem_tags[j]
                node_tags_for_elem = elem_node_tags[j]
                
                # 创建多边形
                poly = vtk.vtkPolygon()
                poly.GetPointIds().SetNumberOfIds(len(node_tags_for_elem))
                
                for k, node_tag in enumerate(node_tags_for_elem):
                    if node_tag in node_tag_to_index:
                        poly.GetPointIds().SetId(k, node_tag_to_index[node_tag])
                
                polys.InsertNextCell(poly)
                
                # 为单元分配FaceID（简化处理）
                face_id_array.InsertNextValue(0)
                cell_count += 1
    
    # 创建PolyData
    polydata = vtk.vtkPolyData()
    polydata.SetPoints(points)
    polydata.SetPolys(polys)
    polydata.GetCellData().AddArray(face_id_array)
    
    # 保存VTK文件
    writer = vtk.vtkPolyDataWriter()
    writer.SetFileName(gmsh_vtk_file)
    writer.SetInputData(polydata)
    writer.Write()
    
    print(f"GMSH VTK模型已保存到: {gmsh_vtk_file}")
    return gmsh_vtk_file


def main_gmsh_processing_with_feature_matching(iges_file_path, faceid_to_features_file, faceid_to_meshsize_file, output_dir="output"):
    """
    GMSH主处理流程 - 使用faceid-to-features数据进行几何特征匹配
    """
    import os
    
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
    
    # 加载数据
    faceid_to_features, faceid_to_meshsize = load_feature_and_meshsize_data(faceid_to_features_file, faceid_to_meshsize_file)
    
    # 初始化GMSH环境
    setup_gmsh_environment()
    
    try:
        # 执行几何特征匹配
        perform_feature_matching(iges_file_path, faceid_to_features, faceid_to_meshsize)
        
        # 配置网格参数
        configure_mesh_parameters()
        
        # 生成并保存网格
        mesh_file_path = generate_and_save_mesh(output_dir)
        
        # 转换网格格式
        try:
            result_file = convert_mesh_to_vtu_with_meshio(mesh_file_path, output_dir)
        except ImportError:
            result_file = convert_mesh_to_vtk_fallback(mesh_file_path, output_dir)
        
        return result_file
        
    except Exception as e:
        print(f"GMSH处理失败: {e}")
        import traceback
        traceback.print_exc()
        return None
    
    finally:
        gmsh.finalize()


# 注释掉可视化相关方法，已移至gmsh_visualization模块
# def visualize_gmsh_physical_groups_2d(msh_file_path, title="GMSH物理组可视化", fig=None, ax=None):
#     """
#     使用matplotlib二维显示GMSH msh文件的物理组颜色映射
#     直接读取msh文件，无需转换为VTK格式
#     """
#     import gmsh
#     import matplotlib.pyplot as plt
#     import numpy as np
#     import matplotlib.patches as patches
#     
#     # 初始化Gmsh
#     gmsh.initialize()
#     gmsh.option.setNumber("General.Terminal", 0)  # 关闭终端输出
#     
#     try:
#         # 加载msh文件
#         gmsh.open(msh_file_path)
#         print(f"已加载GMSH网格文件: {msh_file_path}")
#         
#         # 获取所有物理组
#         physical_groups = gmsh.model.getPhysicalGroups()
#         print(f"找到 {len(physical_groups)} 个物理组")
#         
#         # 创建图形（二维显示）
#         if fig is None or ax is None:
#             fig, ax = plt.subplots(figsize=(12, 8))
#         
#         # 为每个物理组分配颜色
#         colors = plt.cm.tab20(np.linspace(0, 1, len(physical_groups)))
#         
#         # 处理每个物理组
#         for i, (dim, tag) in enumerate(physical_groups):
#             # 获取物理组名称
#             name = gmsh.model.getPhysicalName(dim, tag)
#             print(f"物理组 {i}: 维度={dim}, 标签={tag}, 名称='{name}'")
#             
#             # 获取物理组包含的实体
#             entities = gmsh.model.getEntitiesForPhysicalGroup(dim, tag)
#             
#             # 处理每个实体
#             for entity in entities:
#                 # 获取实体上的单元
#                 element_types, element_tags, element_node_tags = gmsh.model.mesh.getElements(dim, entity)
#                 
#                 # 处理每种单元类型
#                 for elem_type, elem_tags, elem_node_tags in zip(element_types, element_tags, element_node_tags):
#                     # 只处理面单元（三角形和四边形）
#                     if elem_type in [2, 3]:  # 2=三角形, 3=四边形
#                         print(f"处理单元类型 {elem_type}, 单元数量: {len(elem_tags)}")
#                         
#                         # 根据单元类型确定每个单元的节点数
#                         if elem_type == 2:  # 三角形，3个节点
#                             nodes_per_element = 3
#                         elif elem_type == 3:  # 四边形，4个节点
#                             nodes_per_element = 4
#                         
#                         # 处理每个单元
#                         for j in range(len(elem_tags)):
#                             # 计算节点标签在elem_node_tags数组中的起始位置
#                             start_idx = j * nodes_per_element
#                             end_idx = start_idx + nodes_per_element
#                             
#                             if end_idx <= len(elem_node_tags):
#                                 # 获取该单元的节点标签
#                                 node_tags = elem_node_tags[start_idx:end_idx]
#                                 
#                                 # 获取节点坐标
#                                 node_coords = []
#                                 for node_tag in node_tags:
#                                     coord_result = gmsh.model.mesh.getNode(node_tag)
#                                     if coord_result is not None:
#                                         # coord_result是一个元组，第一个元素是坐标数组
#                                         coord_array = coord_result[0]
#                                         if len(coord_array) == 3:
#                                             node_coords.append(coord_array)
#                                         else:
#                                             # 如果坐标不是三维，填充为三维
#                                             if len(coord_array) == 2:
#                                                 node_coords.append([coord_array[0], coord_array[1], 0.0])
#                                             else:
#                                                 node_coords.append([0.0, 0.0, 0.0])
#                                 
#                                 if len(node_coords) == 0:
#                                     continue
#                                     
#                                 # 转换为numpy数组，确保所有坐标都是三维的
#                                 node_coords = np.array(node_coords)
#                                 
#                                 # 提取XY坐标（二维显示）
#                                 x = node_coords[:, 0]
#                                 y = node_coords[:, 1]
#                                 
#                                 # 绘制单元
#                                 if elem_type == 2:  # 三角形
#                                     triangle = patches.Polygon(np.column_stack([x, y]), 
#                                                              closed=True, 
#                                                              fill=True, 
#                                                              color=colors[i], 
#                                                              alpha=0.7,
#                                                              edgecolor='black',
#                                                              linewidth=0.5)
#                                     ax.add_patch(triangle)
#                                 
#                                 elif elem_type == 3:  # 四边形
#                                     quad = patches.Polygon(np.column_stack([x, y]), 
#                                                           closed=True, 
#                                                           fill=True, 
#                                                           color=colors[i], 
#                                                           alpha=0.7,
#                                                           edgecolor='black',
#                                                           linewidth=0.5)
#                                     ax.add_patch(quad)
#         
#         # 设置图形属性
#         ax.set_title(title, fontsize=14)
#         ax.set_xlabel('X坐标')
#         ax.set_ylabel('Y坐标')
#         ax.grid(True, alpha=0.3)
#         ax.set_aspect('equal')
#         
#         # 创建图例
#         legend_patches = []
#         for i, (dim, tag) in enumerate(physical_groups):
#             name = gmsh.model.getPhysicalName(dim, tag)
#             patch = patches.Patch(color=colors[i], label=f'{name} (dim={dim}, tag={tag})')
#             legend_patches.append(patch)
#         
#         ax.legend(handles=legend_patches, loc='upper right', bbox_to_anchor=(1.15, 1))
#         
#         # 自动调整坐标轴范围
#         ax.autoscale_view()
#         
#         print(f"{title} 已使用matplotlib二维显示")
#         
#         return fig, ax
#         
#     except Exception as e:
#         print(f"GMSH二维可视化失败: {e}")
#         import traceback
#         traceback.print_exc()
#         return None, None
#     
#     finally:
#         # 清理Gmsh
#         gmsh.finalize()
# 
# 
# def compare_gmsh_vtk_visualization(msh_file_path, vtk_file_path=None):
#     """
#     同时显示GMSH物理组可视化和VTK FaceID可视化进行对比
#     """
#     import matplotlib.pyplot as plt
#     import os
#     
#     # 创建子图布局
#     fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
#     
#     print("=" * 60)
#     print("开始对比可视化")
#     print("=" * 60)
#     
#     # 显示GMSH物理组可视化
#     if os.path.exists(msh_file_path):
#         print(f"✓ 加载GMSH文件: {msh_file_path}")
#         fig_gmsh, ax_gmsh = visualize_gmsh_physical_groups_2d(msh_file_path, "GMSH物理组可视化", fig, ax1)
#         if fig_gmsh is not None:
#             print("✓ GMSH可视化完成")
#         else:
#             print("✗ GMSH可视化失败")
#             ax1.text(0.5, 0.5, 'GMSH可视化失败', ha='center', va='center', transform=ax1.transAxes)
#     else:
#         print(f"✗ GMSH文件不存在: {msh_file_path}")
#         ax1.text(0.5, 0.5, 'GMSH文件不存在', ha='center', va='center', transform=ax1.transAxes)
#     
#     # 显示VTK FaceID可视化
#     if vtk_file_path and os.path.exists(vtk_file_path):
#         print(f"✓ 加载VTK文件: {vtk_file_path}")
#         fig_vtk, ax_vtk = visualize_vtk_face_ids_matplotlib_2d(vtk_file_path, "VTK FaceID可视化", fig, ax2)
#         if fig_vtk is not None:
#             print("✓ VTK可视化完成")
#         else:
#             print("✗ VTK可视化失败")
#             ax2.text(0.5, 0.5, 'VTK可视化失败', ha='center', va='center', transform=ax2.transAxes)
#     else:
#         print("✗ VTK文件未提供或不存在")
#         ax2.text(0.5, 0.5, 'VTK文件未提供', ha='center', va='center', transform=ax2.transAxes)
#     
#     # 调整布局
#     plt.tight_layout()
#     
#     print("=" * 60)
#     print("对比可视化完成，正在显示图形...")
#     print("=" * 60)
#     
#     # 显示图形
#     plt.show()
#     
#     return fig, (ax1, ax2)